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SUMMARY 

This paper deals with the spectral element modeling of seismic wave propagation at the 
global scale. Two aspects relevant to low-frequency studies are particularly emphasized. 
First, the method is generalized beyond the Cowling approximation in order to fully ac- 
count for the effects of self-gravitation. In particular, the perturbation of the gravity field 
outside the Earth is handled by a projection of the spectral element solution onto the ba- 
sis of spherical harmonics. Second, we propose a new formulation inside the fluid which 
allows to account for an arbitrary density stratification. It is based upon a decomposition 
of the displacement into two scalar potentials, and results in a fully explicit fluid-solid cou- 
pling strategy. The implementation of the method is carefully detailed and its accuracy is 
demonstrated through a series of benchmark tests. 

Key words: Brunt-Vaisala frequency - elastodynamics ~ global seismology - numerical 
modeling - self-gravitation - spectral element method - synthetic seismograms. 
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1 INTRODUCTION 



It has bee n recently estab. 




2002 



Capdeville et al 



i shed b y several auth o rs (|Chaliubl pOOOT ) 



{2m) 



Chaliiib et al 



Komatitsch fc Tromp 



( 20031 )) that the spectral element method 



(SEM) provides an efficient solution to the issue of computing synthetic seismograms in three 
dimensional (3D) models of the Earth. Whereas most of current spectral element studies aim 
at pushing calculations toward high frequencies, where the methods traditionally used at the 
global scale reach their limits, this paper focuses on some physical effects that are critical for 
the lower part of the seismic frequency band: (i) the full treatment of self-gravitation and (ii) 
the ability to take into account any density stratification in the fluid regions of the Earth. 

The ffist novelty of this paper stands in the incorporation of self-gravitation, the effect of 
which is important for seismic and gravimetric observations with periods larger than 100 s. 
All the previously mentioned studie s based upon t he SEM accounted for the effects of gravity 



within the Cowling approximation (jCowling 



194lh . i.e. by neglecting the perturbation of the 



gravity field by seismic waves. The main reason for making this assumption lies in the intrinsic 
difficulty of the problem. Considering the full effects of self-gravitation requires, indeed, to solve 
Poisson's equation for the perturbed gravitational potential which is defined over the whole 
space. Unlike spherical harmonics approaches, the use of a grid-based method such as the SEM 
does not provide a natural framework for the resolution of the exterior problem. Grid-based 
approximations in unbounded domains proceed ffist by restricting the computational domain, 
then by imposing an appropriate condition on the truncating boundary. Different methods 
arise depending on whether the artificial boundary condition (ABC) is local or not. Methods 
based upon a local ABC have the advantage of being computationally inexpensive and valid for 



arbitr a ry geometries. An example o f such methods is the infinite element method {e.g. 



(11992), 



Bf^ttess 



Gerdes fc Demkowica (|l99a )). in which the behaviour of the exterior solution is enforced 
in the radial direction. The second class of methods, based upon a non-local ABC, are not as 
general since they usually require the knowledge of an analytical, or semi-analytical, solution 
to the exterior problem. As a consequence they have very attractive properties regarding their 
accuracy while being restricted to simple (usually spherical) geometries. The non-local ABC can 
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be implemented into the finite el ement metho d within the rigorous framework of a Dirichlet-to- 



Neumann {DtN) operator (e.g. iGivoli 



( 1992I )). This is the approach that we retain here. The 



DtN operator that suits our problem relies on the spherical harmonic deco mposition of the 



solutio n of Laplace's equation outside the Earth. Unlike the one introduced by 



Capdeville et al. 



f)2003^ to couple a time-dependent spectral element calculation to a modal solution in the 
frequency domain, our DtN operator is much simpler to derive because it is applied to a static 
problem. The spectral element discretization of the Poisson-Laplace equation yields a symmetric 
algebraic system which has to be inverted at each time step to obtain the perturbation of the 
gravitational potential. In practice, this is done by iterating a conjugate gradient method, the 
preconditioning of which is critical to carry out routine calculations. 

The other aspect we consider in great detail is the treatment of the fluid part of the Earth's 
core. A parameter which is of particular importance with regard to core dynamics is the squared 
Brunt-Vaisala frequency A^^ that characterizes the local response of the fluid to perturbations 
in density. To first order, the core can be considered as neutrally stratified, i.e. N"^ = 0, because 
a neutral buoyancy is expected in the bulk of a region subject to vigorous convection. However, 
there is seismological evidence for a negative N"^ at the top o f the cor e and a positive N"^ at its 



bottom, with absolute values that can reach 10"^ rad^ ■ s'^ (jMasterj ()l979[ ). Valette & Lesage 
(unpublished)). For the sake of generality, our description of the core's structure will make no 
assumption on the profile of the buoyancy frequency. To this end, we introduce a two-potential 
fo rmulation of the wave equatioii in th e fluid that gen e ralizes the neutral buoyancy formulation 



of 



Komatitsch fc TrompI (l2002br i and 



Chaliub et al 



( 2003^ . Contrary to these studies, that 



considered the velocity potential in the fluid, our decomposition is applied to the displacement 
field in order to obtain natural solid-fiuid boundary conditions for the perturbed gravitational 
potential. An attractive consequence of this choice is to yield a fully explicit sohd-fluid coupling 



strategy, as opposed to the studies mentioned above. Note final 



to the two-potential description proposed by 



Wu Rochester 



y that our formulation is close 



( 19901 ) in the context of core 



dynamics studies, which is optimal with respect to the number of unknowns in the fiuid regions. 
The remainder of the paper is organized as follows. In section |21 we recall the equations of 
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motion in a self-gravitating Earth in their strong and weak form, successively. We introduce 
in particular the two-potential decomposition of the displacement field in the fluid regions and 
we define the DtN operator that permits to handle the equations within a finite domain. In 
section El we recall the principles of the spectral element approximation in space and we make 
a detailed presentation of our explicit time marching algorithm. Finally, numerical results are 
shown in section 0] for a set of spherically symmetric models that validate the implementation 
of the method. 

2 WAVE EQUATION IN A SELF-GRAVITATING EARTH 

In this section we recall the strong and weak forms of the wave equation, which is obtained 
through a first order Lagrangian perturbation around a non-rotating, hydrostatically pre- 
stressed, state of equilibrium. Throughout the paper, the Earth is denoted by © and its outer 
boundary by d^. The solid (resp. fluid) parts of © are referred to as ©5 (resp. ©f), and the set 
of all solid-fluid interfaces is denoted T^sf- Whenever topography or ellipticity is considered on 
9®, B will denote a ball of radius b that contains the aspherical Earth {i.e. , © C i3) and S will 
stand for its spherical boundary (5 = dB). 

2.1 Strong form 

Solving the wave equation within the previous assumptions consists in finding the Lagrangian 
perturbation of the displacement, u, such that: 

ii + ^(u) = -f , (1) 

P 

pAin) = -V-T(u) - V(pu-g) + {V-(pu)}g + pV^, (2) 

where A is the elastic-gravitational operator, T(u) is the Lagrangian incremental stress tensor, p 
is density, g is the acceleration due to gravity, ip is the Eulerian perturbation of the gravitational 
potential, also known as the mass redistribution potential (MRP), and f is the forcing term. As 
usual, a dot over a symbol implies time derivation and Vr (resp. V ■ r) stands for the gradient 
(resp. the divergence) of a given tensor field r. 
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In the (inviscid) fluid regions the stress tensor takes the form: 



T(u) 



pc^V-u I 



(3) 



where c is the speed of sound and I denotes the second-order identity tensor. Neglecting any 
source term in the fluid, the wave equation can then be rewritten as: 



u 



+ V-u s , 



(4) 



(5) 



where s is defined by: 

_ Vp g 

p 

and can be shown to be proportional to the gradient of specific entropy. Another parameter of 
interest in the fluid is the square of the Brunt- Vaisala frequency A^^, which is related to s by: 

A2 = s-g = ^ (Vp - ^g) -g . (6) 

The Brunt- Vaisala frequency arises naturally when analyzing the local stability of the flui d 



since it provides a simple way to formulate the Schwarzschild criterion (jSchwarzschilc 



mm 



An inspection of the expression of the energy reve als, indeed, that the l ocal c oi ivective 



Friedma,n Schut?! ()l978l) 



stabili ty 



Valettel (|l98(tl). 



of the fluid is determined by the sign of A^^ {^-g- 
Actually, A^^ controls the non-seismic part of the spectrum of the elastic-gravitational operator, 

CTeiA): 

aM) = [Min(0,A,L),Max(0,Ar2j] , (7) 




. This implies that 



where A^^^f and N^^^ stand for the extrema of A^ over (Bf 
the corresponding squared eigenfrequencies range in the latter interval. In the Earth, these 
eigenfrequencies merely exceed 50 /iHz, a value which is well below that of the gravest seismic 
oscillation 0*52. 

In this paper, we only intend to compute the seismic part of the fluid outer core's response, 
which is also affected by the variations of A^^. Taking into account a fluid region within the 
framework of the flnite element method is known to be a difficult problem, due to the possi- 
ble splitting of the zero eigenfrequency induced by the numerical discretization of the elastic 
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operator (jHamdi et al, 



19781 ). A key issue to produce a numerical solution free of spurious 



modes is the correct representa tion of the null space of the elastic-gravitational operator, //{A) 



(Rermudez fc Rodrfgue? 



1994( ). An alternative to the discretization of //{A) is to solve the 
wave equation in the range of the operator, 71{A). To proceed, we note from eq. (H)) that an 
acceptable form for any displacement field in TZ{A) is: 

u = Vx + es, (8) 

where x ^ denote two arbitrary scalar fields. Differentiating twice in time and identifying 
each term with the right-hand-side of eq. (jH), we obtain two scalar wave equations, one for each 
potential: 

? + Vx - g + N^^ 



X 



^ ■ 



(9) 
(10) 



Eventually, the MRP -0 appearing in eqs. Q and (fTTH) is obtained by solving the Poisson-Laplace 
equation over the entire space. This writes: 

-47rG'V-(pu) in ©s , 

-47rGV-(pVx + p^s) in©^, (H) 

outside © , 
where G is the gravitational constant. 



2.2 Boundary conditions 



T he complete set of 



m 



boundary conditions for displacement, traction and MRP can be found 



Dahlen fc TrompI (|1998L p. 104). Here we recall these boundary conditions that concern the 
MRP or involve a solid-fluid interface. 

Let S be a given interface in the medium. The condition that the MRP must be continuous 
across S reads: 

Ms = , (12) 
where [ stands for the jump operator across S, defined in accordance with the unit normal 
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vector n: [i/j]-^ = — ip~ and fi points from the — to the + side. The normal derivative of i/j 
can have a jump which is controlled by: 

[Vij-h]^ = -AnG [pu-n]^ . (13) 

The condition that both traction and normal displacement must be continuous across the solid- 
fluid boundaries writes set of equalities on T^sf'- 

u-n = (Vx + es)-n, (14) 
T(u).n = p'ih. (15) 

Note that to obtain eq. (fTHjl we have used eqs. (jHl), (jHl) and 
2.3 Weak form 

The weak form of the wave equation in the solid regions is obtained after multiplying each side 
of eq. (0) with an admissible displacement field w, then integrating over ©5. This writes: 

(u + ^(u) ; pw)^^ = (f;w)0^ (16) 

where ( ; )0^ stands for the scalar product on ©5-. For example, integrating by parts the 
divergence of the stress tensor in eq. Q yields: 

- (V ■ T(u) ; w)^^ = I T(u) ■ Vw (iV - I T(u) ■ n ■ w , (17) 

where n stands for the unit vector normal to T^sf pointing away from ©5. Note that the 
condition of free traction at the surface of the Earth ^0 is naturally satisfied in eq. ()17|) as we 
have set the corresponding integral to zero. On the contrary, the continuity of traction ()15p 
across the solid-fluid boundaries has to be enforced. To proceed, we simply replace the traction 
vector in the surface integral of eq. ()17|1 with its fluid counterpart: 

-(V-T(u); w)^^ = I T(u)-Vwdy- I pin-wdS. (18) 

The weak form of the wave equation in the fluid regions is obtained similarly after dotting each 
side of eqs. (0) and (fTUI) with admissible potentials C, and x, integrating (possibly by parts) over 
(Bf, then forcing the continuity of the normal displacement p4|l across the fluid-solid interfaces. 
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One gets: 



\UdV=-f iVx + Cs)-VidV+ [ u-hCdS, (19) 



and 



\xxdV= [ \ f^ + Vx-g + N^C-ij) xdV . (20) 



In eq. (I19|) . fi denotes the unit vector normal to that points outward the fluid. Note that 
the scaling factor c"^ has been artificially included in eq. (j2n|l in order to get the same left 
hand side as in eq. (fT^ . This will make the description of the time marching algorithm easier 
in section ini 

Now, in order to establish the weak form of eq. (fTT|) . it is convenient to first consider Poisson's 
equation within the finite (spherical) volume B. Multiplying with an admissible potential i/j 
defined over B, then integrating by parts the Laplacian and the divergence we get: 



Vip-VipdV- V^lj-nipdS= (21) 

B Js 

-AttgI f pu-Vi) dV - f pu-ni) dS + f p (Vx + ^ s) ■ rf^ | , 

L>'©S "''5 >'©F J 

with the boundary term involving the normal displacement being null, except in the absence of 
topography {i.e. when S = 9^). It is important to note that the jump condition (fT^ across the 
solid-fluid interfaces is naturally taken into account in (jUJ). This property, which stems from 
the potential decomposition (jSI), is a key argument that guided our choice to work with the 
displacement field (and not the velocity) in the fluid. 

2.4 DtN operator 

The harmonic behaviour of ip outside B has not been considered yet. In order to proceed, let ifj^^^ 
denote the MRP interior to B. At the (spherical) surface cS , consider the expansion of V'mt onto 



the orthonormal basis of real spherical harmonics (see 

oo / 



Dahlen fc TrompI (|l998t ). p.851): 



^Ub,0,^) = E E i^LT{b)yim{e,^) , (22) 
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where {0,(f) are the spherical coordinates and where i^l'^ib) = /^^i„t3^/m dS. It is straightfor- 
ward to extend i/j^^t continuously to a potential i/j^^^ that satisfies Laplace's equation outside B 
and vanishes at infinity: 



QO I / ?)\ ' + ^ 

1=0 m=-l ^ ^ 



CO I / l + l 

yim{a, , r>b. (23) 



The normal derivative of tp^^t on S is readily obtained by differentiating the previous expression 
with respect to r : 

^ oo I 

V^...- n {b, e, if) = -- + 1) E (^) ^'-(^' ^) • (24) 

1=0 m=-l 

Eq. which relates the normal derivative of the potential to the potential itself is called a 
Dirichlet-to-Neumann (DtN) operator on the spherical boundary S. Its action, which is non- 
local, is rather simple to express in the spherical harmonics basis: it consists in multiplying each 
coefficient with ^y^. Recall that the condition that the normal derivative of a given field is 
proportional to the field at the surface is referred to as a Robin boundary condition. Applying 
the DtN operator is therefore equivalent to imposing a Robin boundary condition on every 
component of the spherical harmonics expansion of the original potential, and this yields a 
well-posed problem. 

Taking into account the jump condition (fT!?|) across S, we can write the final weak form of 
the Poisson-Laplace equation as: 

- [ V^-V^! dV+ [ Vtp^^^-htp dS = (25) 
Jb Js 

A-nclj pvi-V^dV+l p (Vx + ^s) ■ (iy 
with: 

« ^ CO I 

/ v^„.,-n^rf5 = --Y^{i + i) v^/:r(6)^""(&) . (26) 

•^^ 1=0 m=-l 

In practice, the infinite sum present in eq. (j26j) will be limited to angular orders / < l^^^ . 
Note that the effect of the truncation is to apply a Neumann boundary condition to the high 
wavenumber content of the MRP, which according to eq. is asymptotically consistent with 
the behaviour of the MRP outside B. 
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3 NUMERICAL APPROXIMATION 

This section deals with the numerical approximation of the wave equation in a self-gravitating 
Earth, which we achieve in two steps. First, the SEM is applied to the weak form of the equations 
in the space domain. Then a finite difference scheme is used to advance the system in time. 
For the sake of conciseness, details of the method are avoided as much as possible unless this 



preve nts the paper from being self-contain ed. The reader is referred to (iKomatitsch fc Vilotte 
19981 ) and to (|Komatitsch fc Tromptll999f) for a general descripti on of the SEM applied to the 



elastic wave equation, and to 



Komatitsch fc Tromp 



2002ai 



2002^) and (jChaliub et al, 



20031) for 



its extension to global seismology, including its parallel implementation on modern computers 
with distributed memory. 



3.1 Spatial discretization 



3.1.1 Hexahedml Mesh 



The first discretization step consists in decomposing the spheric al Earth into a colle ction of 



20021), where 



non-overlapping hexahedral elements. This process is detailed in (jChaljub et al. 
non- conforming interfaces are introduced to avoid an artificial refinement of the grid with depth. 
Such a strategy allows the refinement (or coarsening) of the mesh to be spatially localized, the 
complexity being related to the continuity requirements between elements that do not match 
across the interfaces. For the sake of simplicity, this paper is restricted to the case of a spherical, 
geometrically conforming mesh such as the one represented in fig. IHl Note that taking into 
account the elliptical figure of the Earth or accounting for surface topography would require in 



the self-gravitating case to extend the mesh outward the artificial boundary S. 



3.1.2 Spectral element method 

Based upon the 3D tiling of the sphere, the MRP {ip) as well as the displacement in the solid 
(u) and the potentials in the fluid (x and ^) are approximated using continuous tensorized 
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polynomials. Note that the continuity of the normal displacement within the fluid regions is 
naturally satisfled in the weak forms ()19|) and ((201) • 

The basis of polynomials used on each spectral element are deflned as the shape functions 
of the collocation points. One of the particularity of the SEM is that the collocation points 
are the so-called Gauss-Lobatto-Legendre points, i.e. the exact same points that are used to 
evaluate the integrals present in the weak form of the equations. One consequence of this choice 
is that the matrix representation of th e scalar product is diagona l, a property that allows to 



design explicit time schemes (see e.g. 



(12221)) 



Komatitsch fc Vilottel (|l998h and 



Komatitsch &: Tromp 



3.2 Time evolution 

The different steps of the spatial discretization yield a system of ordinary differential equations 
in time, which writes: 

M5d(t)+ Ksd{t) + Gi^{t)+CsFm = F(t) (27) 

MF^{t)+ KF{tx){t) +CFsd{t) = (28) 

X{t) + Bf t X, V') it) = (29) 

P^{t) = D (d,|,x)(t) (30) 

In the previous equations, d stands for the displacement vector in the solid regions, F is the 
approximation of the source term and i/', x, $, respectively denote the nodal values of the MRP 
and of the displacement potentials in the fluid. M5 is the mass matrix in the solid regions, 
i.e. the matrix representation of the scalar product weighted by density. Similarly, M^? is 
the matrix representation of the scalar product in the fluid regions weighted by the quantity 
c~^. As outlined before, both matrices are diagonal. and are the stiffness matrices which 
arise from the approximation of the volume integrals in eqs. ()18|) and p9p . The discretization 
of the surface integrals in the latter equations yields the solid-fluid coupling matrices Csf and 
Cps- Bi? arises from the discretization of the right hand side of eq. (j20|l and only involves a 
pointwise operation on ^, ^, Vx and i/j. Finally, G, D and P are the matrix representations of 
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the gradient, divergence and Poisson-Laplace operator, respectively. Note that D contains the 
factor AttGp and that P is symmetric according to eqs. and (j^Uj) . 

To advance the equations for ward in time we use the exphcit, second-order accurate. New- 



mark scheme {e.g. 



Hugues 



19871 ). Let for example X„ denote the snapshot at time t„ of one of 
the unknown vectors d, % or ^ involved in eqs. ()27H29|) . The values of X and its time derivative 



at the next time step are extrapolated as follows: 

Xn+l 



X„ + At X„ H — — X^ 
At / •■ 

Xn H — — ^X„ + X^+i 



(31) 
(32) 



As it is readily seen from the previous equations, the algorithm is fully explicit in terms of X 
and consists in a simple centered finite difference scheme in X. The process of updating the 
time derivatives of X is achieved in two steps: first X„+i is computed from the discrete version 
of the wave equation ()27fEn|l by inverting a diagonal mass matrix (M^ or Mi?), then X„+i can 
be updated using (jH^ . Note that the wave equation has to be solved in the fluid regions first, 
since the coupling operator C^i;' in eq. ^I7\i acts on ^n+i which is not known at time 

Let us stress that the coupling between the fluid and the solid regions does not require 
iterati ons of eqs. (I31I32I) as th i s would be the case if a velocity potential formulation was used 



.e.g. 



Komatitsch et al 



200G; 



Chaliub et al 



20031 ). This attractive property stems from the 



potential decomposition (jHl) applied to the displacement which is the explicit variable in the 
Newmark scheme. 

The previous remark remains valid when the full effects of self-gravitation are taken into 
account. The computation of the MRP from the displacement field is indeed explicit in the sense 
that it does not involve any time derivative X or X. Needless to say, thi s task is expensive a s it 



Deville et al. 



20021) . In 



requires to formally invert the symmetric, ill-conditioned matrix P {e.g. 
practice, we solve eq. for the MRP with a conjugate gradient (CG) method which iterations 
are stopped when the residual is decreased by a factor e to be chosen. The issue of building an 
efficient preconditioner for the Poisson-Laplace solver is not addressed in this paper, but it is 
certainly critical in order to avoid a performance bottleneck. 



Modeling wave propagation in a self- gravitating Earth 13 



4 NUMERICAL RESULTS 

In this section, we demonstrate the vahdity of our approacli tlirough a couple of examples for 
which a reference, semi-analytical, solution can be derived. First, the two potentials formulation 
is tested within the Cowling approximation, i.e. without computing the MRP, for models having 
a constant Brunt- Vaisala frequency. T hen, the effects of mass redistri bution are included in a 



simplified version of the PREM model (jPziewonski fc Anderson 



12m- 



4.1 Validation of the two- potentials formulation 

In order to define some benchmarks to test our formulation, we consider the radial Earth 
model of fig. ^ The model is adapted from PREM, with a smaller number of regions (6 instead 
of 13). In particular, the details of the crustal structure as well as the presence of a global 
ocean are ignored to ease the computation. This reference model is further constrained to fit 
a given profile of the squared Brunt- Vaisala frequency in the fluid outer core. To proceed, we 
simply vary the P- velocity in eq. keeping the density, its gradient and the gravitational 
acceleration unc hanged. Note that a rea listic way would be to adjust density rather than P- 



velocity (see e.g. IWu fc Rochesterl ()l993r )) because the latter is much better constrained in the 
Earth. However, acting on the P-velocity profile is straightforward and still fully acceptable for 
numerical validation purposes. 

Fig.Elshows three models that were built following the above procedure. The 'N' label refers 
to a neutrally stratified outer core {i.e. with A^^ = 0), whereas the models labelled 'S' and 'U' 
correspond to a stable and unstable stratification, respectively. For the sake of simplicity, we 
chose the value of the squared Brunt- Vaisala frequency to be constant throughout models 'S' 
and 'U', respectively equal to A^ = 10^^ rad^ ■ s~^ and A^ = —5 10^^ rad^ ■ s~^. These values 
correspond to the extrem a that are expected from the inversion of seismic free oscillations of 



the Earth ( Master 



(19790, Valette & Lesage, unpublished). Note that the values of A^ within 
PREM are about one order of magnitude smaller, as illustrated by the similarity of the PREM 
P-velocities to those of a neutrally stratified profile. 



All three models are excited by a shallow explosive point source which time dependence 
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is a Ricker wavelet {i.e. the second derivative of a Gaussian bell) with dominant frequency 
/o = 1 niHz. The source is located at one grid-point from the Equator, at latitude Og ~ —1.128° 
and depth ds — 61 km, and the receivers sit along the Equator. Fig. El shows the longitudinal 
displacement recorded at an epicentral distance of 90° in the three models. The traces were 
computed within the Cowling approximation using a summation of the eigenmodes of each 
model. The waveform differences illustrate the sensitivity of the seismic waves to the stratifica- 
tion of the fluid core and suggest that models 'U' and 'S' constitute a demanding benchmark 
for the two potentials formulation. In figs. |3] and the spectral element results obtained in 
those two models are compared to the modal solutions for a couple of epicentral distances. The 
two solutions are in very close agreement with the largest relative differences being as small as 
one per mil over the time interval considered. 

The spectral element grid used to carry out the calculations is shown in fig. IHl It consists of 
640 elements in which the polynomial degree varies from 3 to 10 in the radial direction and is 
kept constant, equal to 8, in the tangential direction. The total number of gridpoints is 334,368 
corresponding to a number of points p er wavelength much greater th an 5, which is the empirical 



ratio to get an accurate solution {e.g. iKomatitsch fc Vilottd ()l998| )). This explains the perfect 



match between the spectral element calculations and the reference solutions. 
4.2 Validation of the whole formulation 

As a last example, we consider the computation of the elastic-gravitational response of the Earth 
model of fig. H This test presents all the difficulties mentioned in this paper: the stratification of 
the fiuid core is arbitrary and the physical description includes the full effects of self-gravitation. 

The parameters of the simulations are slightly different than above, since the source domi- 
nant frequency is set to a graver value /o = 0.5 mHz, and the source latitude is now 9s — —2.64°. 
The spectral element grid is consequently adapted, and roughly coarsened by a factor of two 
in each direction compared to the one of fig. El 

In order to check that the test is demanding enough with regard to the implementation of 
self-gravitation, we compare in fig. [71 the surface longitudinal displacement recorded with or 
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without including the perturbation of the gravitational potential. Both traces were computed 
by normal modes summation and recorded at an epicentral distance of 90° for about 10 hours. 
The differences in phase and amplitude illustrate that the Cowling approximation is not valid 
in the frequency range of the experiment. 

Finally, the results obtained with the SEM are compared to the reference solution in fig. |H1 
Two cases are considered that correspond to a different accuracy of the spectral element solution 
regarding the CG resolution of the discrete Poisson-Laplace equation (j29|l . In the first case the 
CG iterations are stopped when the residual is decreased by three orders of magnitude, which 
means that e = 10^^. The resulting spectral element solution is clearly not accurate enough 
and contains a secular term that seems to break the conservation of energy at the discrete 
level. To correct this behaviour, we consider a second test where the stopping criterion is fixed 
to e = 10-^ In that case, the calculation is stable upon the time interval considered and the 
accuracy of the spectral element solution is found to be acceptable, its relative difference with 
the reference solution being less than a few per mil. 

In each of the previous cases, the angular order truncation in eq. was set to l^^^ = 20, 
based on the a priori knowledge of the dispersion relation in PREM. The effect of underesti- 
mating the truncation order is to add oscillations to the spectral element solution (not shown 
in this paper). It is interesting to note that the two possible sources of numerical errors (e too 
big or /jnj^x too small) lead to a different signature. This provides two different diagnostics that 
permit to build a spectral element solution with arbitrary accuracy. 

5 CONCLUSIONS 

We have shown how the SEM should be adapted to account for two effects relevant to global 
seismology: the full treatment of self-gravitation and the ability to consider any density strat- 
ification in the fluid outer core. The accuracy of the method has been illustrated through a 
series of numerical tests conducted in spherically symmetric models. With the incorporation of 
the two aforementioned effects, we believe the SEM will provide new estimates of the elastic- 
gravitational response of 3D models of the Earth. 
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Figure 1. Variation with depth of density (dashed curve), P- velocity (solid curve) and S- velocity 
(dot-dashed curve) within the E arth-like model used in this paper. The model is adapted from PREM 



( Dziewonski fc Anderson 



198ll ) with the complexity of the lithospheric structure being removed to 



simplify computation. 
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Figure 2. Different profiles of P-velocity used to test tlie two-potentials formulation in the fluid 
outer core. The dashed curve represents the variation of the sound speed within the model detailed 
in fig. Each solid curve corresponds to a modification of that profile such that the square of the 
Brunt- Vaisala frequency is constant throughout the fluid. The label 'N' corresponds to a neutrally 
stratified outer core, whereas 'S' (resp. 'U') stands for a stable (resp. unstable) stratification for which 
Ar2 = io~7 rad^ . s-2 (^esp. N"^ = -5 IQ-*^ rad^ • s'^). 
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Figure 3. Time window of the longitudinal surface displacement recorded at 90° in the models 
labelled 'N', 'S' and 'U' in fig. |2l The large waveform differences stem from the sensitivity of the 
seismic modes to the variation of the P- velocity within the three models. 
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Figure 4. Radial (left panel) and longitudinal (right panel) components of the surface displacement 
recorded at 45° (top) and 90° (bottom) in the model labelled 'S' in fig. [51 In each plot, the spectral 
element solution (dashed line) is compared to the normal modes reference (solid thin line) and the 
residual (solid bold line) is amplified by a factor of 10. 
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Figure 5. Radial (left panel) and longitudinal (right panel) components of the surface displacement 
recorded at 45° (top) and 90° (bottom) in the model labelled 'U' in fig. |21 In each plot, the spectral 
element solution (dashed line) is compared to the normal modes reference (solid thin line) and the 
residual (solid bold line) is amplified by a factor of 10. 
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Figure 6. Spectral element mesh used to compute the results shown in figs. 0] and |S1 Two blocks 
of the 3D mesh have been removed to allow a view inside the volume. The mesh is composed of 640 
spectral elements with varying polynomial ord er, for a total numb er of gridpoints equal to 334,368. 



The process of building the mesh is detailed in (|Chaliub et al 



20031 ). This image was generated using 



the visualization software pV3 (http://raphael.niit.edu/pv3/pv3.html). 
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Figure 7. Longitudinal component of the surface displacement recorded at 90° in the Earth model 
of fig. n The trace computed with the full treatment of self-gravitation (solid thin line) is compared 
to the one computed within the Cowling approximation (dashed bold line). The waveform differ- 
ences illustrate that the effect of the MRP cannot be neglected at the frequencies considered in this 
experiment. 
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Figure 8. Longitudinal surface displacements recorded at 90° in the Earth model of fig. ^ The left 
(resp. right) plot corresponds to a low (resp. high) accuracy test in which the CG iterations used to 
compute the MRP are stopped when the residual is decreased by 3 (resp. 5) orders of magnitude. 
In each plot, the spectral element solution (dashed line) is compared to the normal modes reference 
(solid thin line) and the residual (solid bold line) is amplified by a factor of 10. 



